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Abstract 

The  Lattice  Boltzmann  Method  continues  to  garner  interest  in  fluids 
research,  particularly  with  its  ability  to  accurately  simulate  laminar  flows 
in  the  incompressible  region.  This  interest  can  be  attributed  in  part  to  the 
ease  of  implementation  the  Lattice  Boltzmann  Method  provides;  including 
a  lack  of  complex  differential  terms  and  a  linear  approximation  of  the 
collision  operator  contained  in  the  Boltzmann  equation.  In  this  work,  the 
traditional  Lattice  Boltzmann  solver  is  augmented  with  support  for 
immersed  boundaries,  thermal  flows,  and  microchannel  flows.  Thermal 
and  micro-enabling  support  is  demonstrated  and  validated  through 
Rayleigh  convection  in  a  square  channel  and  thermally  coupled  Poiseuille 
flow  through  a  microchannel,  respectively. 
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1  Introduction 

1.1  Background 

The  use  of  the  Lattice  Boltzmann  Method  (LBM)  for  numerically  simulating 
secondary  (amorphous)  phase  constituents  (most  often  represented  as  a 
laminar,  incompressible  fluid)  within  the  context  of  material  science-based 
applications  is,  for  these  purposes,  strategically  based  on  two  primary 
motivations.  The  first,  and  possibly  the  most  dominant,  relates  to  the 
omnipresent  problem  of  spatiotemporal  scale  compatibility  and  reconcilia¬ 
tion,  wherein  one  attempts  to  make  macroscale  material  decisions  based  on 
microscale  or  nanoscale  performance  criteria.  The  second  motivation 
relates  to  concerns  about  the  accuracy  of  the  method,  particularly  with 
respect  to  flows  over  complex  geometries. 

While  the  details  surrounding  these  two  strategic  purposes,  including  a 
thorough  description  of  the  “pros”  and  “cons”  for  using  LBM  over  the 
more  traditional  Computational  Fluid  Dynamics  (CFD)  techniques  are 
included  within  the  first  report  of  this  series  (Allen  et  al.  2014).  This  first 
report,  in  terms  of  validated  examples  and  algorithm  development, 
provided  only  an  elementary  preview  into  the  vast  capabilities  and 
sophistication  of  LBM.  The  necessary  extensions  leading  ultimately  to  the 
resolution  of  the  assigned  motivating  purposes  were  left  largely  unfulfilled. 
The  goal  therefore,  of  this,  the  second  report  in  the  series,  is  to  fulfill  these 
aforementioned  purposes,  and  advance  current  understanding  through  a 
series  of  verified,  algorithmic  extensions  to  the  method. 

In  Allen  et  al.  2014,  an  atomistic  fluid  solving  method  founded  on  statistical 
mechanics  and  kinetic  theory  was  introduced.  The  LBM  includes  no 
continuum  assumption  with  regard  to  the  flow;  rather  the  fluid  is  described 
by  individual  distribution  functions  that  describe  the  expectation  of  finding 
a  particle  in  some  phase  and  temporal  domain  (Allen  et  al.  2014).  In  this 
model,  macroscopic  quantities  such  as  density  and  velocity  are  easily 
computed  as  velocity  moments  of  the  distribution  functions.  From  these 
moments,  nearly  all  other  quantities  specific  to  the  LBM  can  be  described. 
The  accuracy  of  the  LBM  with  respect  to  macroscopic,  isothermal  fluid 
simulations  against  static,  rigid  boundaries  was  demonstrated  in  Report 
One  through  use  of  a  number  of  benchmark  flow  cases.  It  was  shown  that 
results  of  these  demonstrations  were  found  to  be  in  very  good  agreement 
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with  Navier-Stokes  (NS)  based  solutions  in  all  cases,  ensuring  the  validity  of 
the  LBM  as  an  accurate  fluid  solving  method.  The  Lattice  Boltzmann 
Equation  (LBE)  as  presented  in  Report  One  is: 

ft  (*  +  cA ,  t  +  8t )  ■ -  fi  (x,  t )  =  ^(f.  (x,  t )  -  f;q  (x,  t)) + F&  ( 1) 

However,  as  stated  in  Report  One,  a  determining  factor  for  the 
implementation  of  the  LBM  was  its  ability  to  simulate  microscale  flows,  a 
feature  not  often  available  to  more  traditional  solvers  due  the  to  a 
continuum  assumption  for  the  fluid  model.  Further,  it  was  stated  that  the 
LBM  is  able  to  more  accurately  simulate  flow  against  complex  geometries 
due  the  mesh  being  of  a  simpler  quality  and  less  computationally  expensive. 

As  stated,  this  report  attempts  to  address  these  proposals  by  providing  an 
updated  implementation  of  the  Immersed  Moving  Boundary  (IMB)  method 
made  possible  by  accurately  capturing  the  momentum  exchange  between 
particles  and  fluid  interface.  Also,  this  report  demonstrates  the  capacity  of 
the  method  to  support  microchannel  simulations  by  introducing  an  accurate 
model  of  the  Knudsen  layer  with  Diffuse  Scatter  Boundary  Conditions 
(DSBC),  wherein  slip  velocity  on  the  channel  walls  is  introduced.  With 
regard  to  the  IMB  method,  certain  properties  were  desired;  specifically,  that 
the  flow  is  able  to  respond  to  irregular  geometries  despite  the  lattice  domain 
property  of  the  LBM  that  introduces  stair  step  geometries.  Initial  rigid 
boundary  implementation  should  be  changed  to  a  dynamic  one,  able  to 
deform  per  particle  forcing  interactions.  With  such  an  implementation, 
complex  channels  and  boundaries  can  be  accurately  represented,  as  are 
often  encountered  in  computation  fluid  domains. 

With  regard  to  the  microscale  implementation,  the  chosen  method  should 
accurately  account  for  fluid  properties  due  to  rarefaction  effects  at  the 
microscale.  Among  these  effects,  is  the  presence  of  Knudsen  layer  whose 
existence,  in  many  cases,  dominates  fluid  interactions.  Due  to  higher 
collision  frequency  at  the  wall,  the  implementation  should  also  be  able  to 
simulate  a  wall  slip  velocity  condition.  It  is  to  be  shown  that  the  DSBC, 
whose  roots  lie  in  kinetic  theory,  is  suitable  for  such  conditions.  Further, 
implementation  of  the  tangential  momentum  accommodation  coefficient, 
being  a  property  inherent  to  the  material,  allows  the  DSBC  to  be  extended 
for  modeling  not  only  the  aforementioned  complex  geometries,  but  also 
various  material  types. 
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1.2  Objective 

This  report  seeks  to  demonstrate  that  the  accuracy  of  the  LBM  is  main¬ 
tained  even  with  the  aforementioned,  implemented  extensions.  While 
Navier-Stokes  based  solvers  tend  to  dominate  CFD  related  simulations,  the 
LBM  has  received  much  interest  in  recent  years,  due  in  part  to  being  able  to 
produce  accurate  results  where  Navier-Stokes  approximations  fail; 
particularly  as  a  the  modeled  scale  decreases  (Allen  et  al.  2014).  In-house 
developed  code  compared  against  international  research  efforts  demon¬ 
strate  the  validity  of  the  current  LBM  implementation  through  a  variety  of 
popular  simulation  cases,  including;  drag  coefficient  calculations  for  an 
immersed  cylinder,  Rayleigh  convection  in  a  square  channel,  and  microscale 
Poiseuille  flow. 

1.3  Approach 

This  second  report  is  intended  to  serve  as  an  extension  to  Report  One.  It 
was  shown  in  the  first  report  (Allen  et  al.  2014)  that  the  LBM  can  return 
the  Navier-Stokes  equations  using  the  Chapman-Enskog  expansion.  While 
this  report  does  not  return  implemented  contributions  to  the  NS 
equations,  certain  NS  contributions  can  be  seen  in  the  thermal  support 
equations.  Further,  microscale  implementations  step  outside  the  realm  of 
the  macroscopically  employed  NS  scenarios  and  employ  concepts  derived 
from  kinetic  theory,  being  central  to  the  foundations  of  the  LBM. 
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2  Immersed  Moving  Boundary 

In  this  work,  support  for  the  IMB  is  implemented  in  the  LBM  to  simulate 
particulate  flows  for  complex  moving  geometries.  The  IMB  was  originally 
developed  by  Peskin  (2002)  for  biological  simulation  of  blood  flow  around 
the  heart  and  has  since  quickly  evolved  into  an  effective  method  for  fluid¬ 
boundary  simulation  for  incompressible  flows.  This  section  is  validated  by 
applying  laminar  flow  around  an  immersed  cylinder.  Like  the 
implementation  of  the  classical  LBM,  the  IMB  method  follows  a  similar  ease 
of  implementation.  No  complex  or  computationally  intensive  differential 
terms  are  needed  and  the  translation  from  Lagrangian  space  to  Eulerian 
space  is  relatively  straightforward.  Further,  implementation  results  of  the 
IMB  with  the  single  relaxation  time  Bhatnagar-Gross-Krook  (BGK)  collision 
remain  comparable  to  multi-relaxation  models  (Niu  et  al.  2006). 

In  the  IMB  approach,  boundary  nodes  are  represented  as  a  mesh  of 
underlying  Lagrangian  markers  able  to  be  moved  in  continuous  space. 

This  contrasts  with  the  discretized  Eulerian  markers  that  mark  the  flow  as 
shown  in  Figure  1.  Once  the  flow  on  the  Eulerian  markers  reaches  a 
Lagrangian  boundary  marker,  a  calculated  force  is  exerted  on  the 
Lagrangian  boundary  marker,  accelerating  the  boundary  and  updating  its 
translational  and  angular  velocities.  The  exerted  force  on  the  Lagrangian 
markers  is  then  pushed  back  on  the  Eulerian  markers  using  a  carefully 
chosen  mollifier  to  ensure  agreement  with  Newton's  Third  Law  (Feng  and 
Michaelides  2003;  Peskin  2002).  Functions  that  maybe  categorized  as 
mollifiers  must  be  smooth  and  the  functionally  bounded  area  must  be 
unity.  This  makes  a  well-chosen  mollifier  ideal  for  reapplying  the 
boundary  force  onto  the  flow. 

The  IMB  algorithm  procedure  occurs  as  follows: 

1.  Collide  the  particles. 

2.  Add  calculated  forces  to  the  collisions. 

3 .  Stream  updated  distribution  functions . 

4.  Update  Lagrange  markers  using  Lagrange  polynomials. 

5 .  Calculate  bounceback  contribution  at  Lagrange  points. 

6.  Update  Eulerian  mesh  using  mollifier. 

7.  Repeat  steps  1-6  until  convergence. 


ERDC  TR-14-6;  Report  2 


5 


Figure  1.  Non-rigid  Lagrangian  mesh  imposed  over  Eulerian  mesh. 
The  dark,  hollow  points  represent  Langrangian  nodes,  while 
vertices  of  coordinate  lines  represent  Eulerian  nodes. 


As  shown  by  Equation,  l,  the  LBM  allows  the  addition  of  external  forces.  In 
the  IMB,  the  external  force  is  the  force  that  is  exerted  by  the  Lagrangian 
mesh  onto  the  Eulerian  mesh.  While  there  are  many  schemes  to  add  an 
external  force  to  the  LBM  equations,  in  this  work,  the  force  adding  scheme 
demonstrated  by  Mohamad  et  al.  (2010),  as  shown  in  Equation  2  is  applied. 


(2) 


This  requires  the  force  to  be  split  into  eight  separate  elements  (F0  is  always 
zero  due  to  c0)  to  be  added  to  the  respective  distribution  functions.  As 
previously  mentioned,  the  single  relaxation  time  BGK  collision  term  is 
used  throughout  this  work.  After  colliding  the  particles,  the  force  is  then 
added  as  previously  shown.  The  particles  are  then  streamed  across  the 
Eulerian  markers. 

To  represent  the  boundary  locations,  Lagrangian  markers  are  placed  at  Xs 
for  s  e  H  where  H  is  the  fluid  domain,  s  represents  the  sth  marker,  and  Xs  is 
the  position  of  that  marker  at  AS(AS,  Es).  Note  that,  in  this  work,  capital 
letter  quantities  are  used  to  represent  Lagrangian  quantities  while 
lowercase  letter  quantities  represent  Eulerian  quantities. 
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To  calculate  the  distribution  functions  at  a  Lagrangian  marker,  two- 
dimensional  Lagrange  interpolation  is  employed  on  the  surrounding 
Eulerian  markers  by  Equation  3.  Here,  interpolate  all  Lagrangian  marker 
distribution  functions  accurate  to  the  third  order. 


•  \ 

hot  X  — r 

rr  xhi 

Ys-yim 

V 

11  v 

xij  xkj 

m=l,m^j  Uij  }Jim 

/afef 


(3) 


Once  the  distribution  functions  have  been  interpolated,  they  are  updated 
using  the  bounceback  scheme  (Zou  and  He  1997)  presented  in  Equation  4 
to  ensure  a  no-slip  condition  on  the  boundary,  where  /?  is  the  opposite 
direction  of  a  or  Cp  —  —  caand  utot  is  the  total  velocity  written  in  terms  of 

the  translational  velocity,  u0,  and  the  angular  velocity,  <u0,  shown  in 
Equation  4.1. 


f,(x„t+s,)=f'(x„t) 


(4) 


where 


utot  =  u„  +  c dnx[X'-X, 


(4.1) 


The  force  density,  FS(XS,  t)  at  Xs  that  is  to  be  exerted  on  the  nearby 
Eulerian  markers  is  calculated  by  Equation  5.  The  difference  in  the 
opposite  distribution  functions  is  representative  of  a  momentum  exchange 
between  the  particles  that  results  in  a  change  in  the  boundary's  velocity. 
The  direction  of  Fs  is  determined  by  the  magnitude  of  this  exchange  when 
multiplied  by  the  lattice  velocity. 


(5) 


With  FS[XS,  t )  calculated,  the  force  can  now  be  distributed  to  the  flow 
through  the  use  of  the  aforementioned  mollifier.  The  total  force 
redistributed  to  the  flow  is  given  by  Equation  6. 


Fr(r,t)  =  jFs(x„t)s(r-X)ds 


(6) 
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where  8h  is  the  chosen  mollifier.  In  his  work,  Peskin  (2002)  presented  two 
delta  functions  for  8h,  shown  in  Equation  7  and  Equation  8. 


.  \  \ 

1  .  JVC  I  I 

.  v  —  1  +  cos  —  ,  x  <2 

«»M=W  2  JJ’1 

0,  otherwise 


Applying  8h(x )  to  both  the  x  and  y  Eulerian  coordinates  yields  Equation  9. 
The  effective  force  distribution  geometry  is  a  square  where  the  Lagrangian 
marker  is  located  at  the  center  and  the  length  of  the  sides  is  two  for  both 
Equation  7  and  Equation  8.  In  this  work,  Equation  8  is  used  exclusively 
primarily  due  to  ease  of  implementation  without  loss  of  accuracy 
(Appendix  A). 


(9) 


Independent  of  the  mollifier,  because  the  fluid  domain  is  discretized  into 
Eulerian  markers  with  constant  spacing  8X,  Equation  6  must  be  discretized 
accordingly.  (It  is  important  to  note  that  summing  the  mollifier  on  discrete 
steps  of  Ax  still  sums  to  unity,  proving  that  no  magnitude  of  force  is  lost  in 
redistribution  to  the  flow.  See  Appendix  A  for  a  more  formal  proof  of  this). 
Placing  Equation  9  into  Equation  6  yields  the  final  force  distribution. 


(10) 
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2.1  Object  properties 

In  this  implementation  of  the  IMB,  Lagrange  markers  are  initialized  with  a 
link  to  the  nearest  Lagrange  markers  for  quick  access  to  information 
contained  by  surrounding  Lagrange  nodes.  Following  this,  it  is  possible  to 
start  at  an  arbitrary  Lagrange  marker  and  follow  its  links  to  many  other 
Lagrange  nodes.  The  term  object  is  applied  to  the  collection  of  such  a 
group  of  linked  Lagrange  markers.  Further,  this  linking  scheme  allows  the 
possibility  to  create  as  many  objects  as  desired  which  are  able  to  interact 
with  other  objects  with  little  description  of  object  interaction. 

In  this  implementation,  each  Lagrange  marker  is  weighted  by  a  mass  ms, 
such  that  an  object's  mass  is :m0  —  £s  ms.  The  first  moment  of  mass 
divided  by  m0  is  shown  in  Equation  n  and  provides  the  location  of  the 
center  of  mass  of  the  object. 


*_=— (in 

m 

"Lo  S 

To  obtain  the  total  force  on  an  object,  simply  sum  the  balance  force  of  the 
respective  Lagrange  markers'  forces  contained  by  the  object  as  shown  in 
Equation  12. 


K  = 

S 


(12) 


This  now  allows  the  opportunity  to  turn  to  Newton's  Second  Law  to 
calculate  the  translational  velocity  of  an  object,  provided  in  Equation  13. 


du,, 


Fdt 


m„ 


(13) 


Similarly,  the  total  torque  on  the  object  is  given  by  Equation  14.  Note  that 
the  torque  acts  in  the  third  dimension,  and  is  not  present  in  the  2- 
dimentional  lattice  (D2)  with  nine  discrete  velocities  (D2Q9)  model  (Allen 
et  al.  2014). 


(i4) 

S 

The  moment  of  inertia,  I0,  is  given  by 
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I0  =  ms(Xs-Xcmf  (15) 

From  the  object’s  torque  and  moment  of  inertia,  the  angular  velocity  can 
be  obtained,  which  is  shown  in  Equation  16. 


(16) 


Finally,  considering  both  the  translational  and  angular  velocity 
contributions,  the  total  velocity  is  calculated  by  Equation  4.1. 

In  this  implementation,  the  object  is  able  to  rotate  about  its  center  of 
mass.  This  is  accomplished  by  employing  the  rotation  matrix  shown  in 
Equation  17,  where  6  is  given  in  Equation  18. 

cos  6  sin0 
-sin0  cos  (9 

d0  =  c6odt  (18) 


The  parameter  X's  denotes  the  new  location  of  Lagrange  marker  Xs  after 
rotation  and  is  given  by  Equation  19. 

Xs=R(Xs-XcJ  +  Xcm  (19) 

2.2  Validation 

To  validate  the  computational  accuracy  of  the  IMB  method,  laminar  flow 
across  an  immersed  cylinder  is  considered.  This  flow  problem  has  been 
studied  extensively  and  has  been  shown  to  preserve  accuracy  in  the  LBM. 
(Niu  et  al.  2006).  However,  this  problem  can  be  difficult  to  model 
effectively  due  to  the  stair  step  geometry  imposed  on  boundary  walls  by 
lattice  discretization.  The  IMB  method  remediates  this  problem  due  to  the 
boundary  being  composed  of  Lagrangian  markers  with  no  ties  to  any 
specified  Eulerian  marker.  One  common  measurement  of  flow  around  an 
immersed  cylinder  is  the  drag  coefficient  CD .  This  quantity  may  be  defined 
in  terms  of  the  streamwise  component  of  the  force  acting  on  the  boundary 

p 

by  CD  =  — where  it0  is  the  inlet  stream  velocity  and  r  is  the  radius  of 

J  PavgU20r’  u  J 

the  cylinder. 
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To  set  up  this  simulation,  the  cylinder  was  initially  placed  at  x  — 

where  L  and  H  are  the  length  and  height  of  the  channel,  respectively.  A 
prescribed  velocity  of  u0  is  set  at  the  inlet.  The  top  and  bottom  walls  are 
treated  by  standard  no-slip  bounceback  conditions  and  the  outlet  density 
is  extrapolated  from  nearby  lattice  sites.  Further,  the  size  of  the  channel  is 
defined  by  L  =  8or  and  H  =  8or.  It  can  be  seen  in  Figure  2  that  two 
vortices  develop  behind  the  object.  The  x-wise  length,  or  recirculation 
length,  of  the  vortices  can  be  expressed  in  terms  of  the  radius  of  the 
cylinder  by  Lw  =  where  Lr  is  the  measured  recirculation  length. 

Results  of  yielded  drag  coefficient  and  recirculation  for  Re  =  20  are 
compared  with  other  LBM  implementations  by  Niu  et  al.  (2006)  and  He  et 
al.  (1998)  and  are  presented  in  Table  1. 


Figure  2.  Flow  around  and  immersed  cylinder  for  Re  =  20. 


Table  1.  Drag  coefficient  and  recirculation  length  comparison  for  Re  =  20. 


C  D 

Abs.  Error 

Lw 

Abs.  Error 

Niu  et  al. 

2.144 

- 

1.89 

- 

He  et  al. 

2.152 

0.373 

1.84 

2.645 

Current 

2.259 

5.363 

1.95 

3.174 
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3  Thermal  Interactions 


Thermal  support  in  the  LBM  has  been  an  area  of  interest  for  accurate 
simulation  of  thermal  particulate  interaction  for  various  scales  and  has  led 
to  the  development  of  thermal  implementation  methods.  Among  these 
methods  are  the  multispeed  method  (McNamara  et  al.  1995),  the  passive 
scalar  method  (Shan  1997)  and  the  thermal  energy  distribution  function 
method  (He  et  al.  1998). 

The  multispeed  method  extends  the  LBM  equation  to  include  calculations 
of  macroscopic  temperature  through  implementation  of  new,  higher  order 
velocity  terms.  While  the  multispeed  approach  may  be  straightforward  for 
the  LBM,  results  from  this  method  are  often  unstable.  Further,  the  range 
of  supported  temperatures  is  limited.  (McNamara  et  al.  1995) 

The  passive  scalar  method  implements  thermal  support  by  introducing  a 
second  distribution  function  to  simulate  thermal  interactions.  This 
method  assumes  that  viscous  heat  dissipation  and  compression  work  done 
by  the  pressure  can  be  ignored,  but  does  suffer  from  the  same  numerical 
instability  as  the  multispeed  method. 

In  this  work,  the  thermal  energy  distribution  function  method  proposed  by 
He  et  al.  (1998)  is  implemented.  Following  the  work  done  by  Peng  et  al. 
(2003),  the  thermal  energy  distribution  function  is  simplified  by  neglecting 
the  compression  work  done  by  the  pressure  and  the  viscous  heat  dissipa¬ 
tion.  Eliminating  these  contributions  is  a  stipulation  also  imposed  by  the 
aforementioned  passive  scalar  method  and  are  justified  in  incompressible 
fluid  interaction  as  described  by  Peng  et  al.  (2003)  Further,  these  contribu¬ 
tions  are  contained  in  a  complex  gradient  term  that  complicates  the  LBM 
and  was  shown  by  Peng  et  al.  (2003)  to  double  computation  time  in  many 
situations. 


The  density  distribution  function  and  thermal  distribution  functions  as  per 
the  work  of  He  et  al.  (1998)  are  presented  in  Equations  20  and  21. 


fi(x  +  ciSt,t  +  St)-fi(x,t)  = 


T„F;6t 


F  +  ^St 


(20) 
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9i  x  +  c&J  +  SA-gAxJ 


+  25t 


gfctygrfat 


Tcfi{x,t)qi8t 

*c  +  |$ 


(21) 


where  /)  and  gt  are  the  new  distribution  functions  introduced  by  He  et  al. 
(1998)  and  are  written  in  terms  of  and  gt,  where  ft  is  the  traditional 
density  distribution  function  and  gt  is  the  thermal  distribution  function  and 
qL  represents  the  viscous  heating  and  compression  work.  The  respective 
relaxation  times  are  now  denoted  as  rv  and  rc  for  the  density  and  thermal 
distribution  functions,  respectively.  Also  proposed  by  He  et  al.  (1998)  are 
related,  but  simpler  functions,  shown  in  Equations  22  and  23. 


where  Ft  is  the  distribution  of  a  force  G  acting  on  the  system  and  the  force 
adding  scheme  is  applied  in  the  collision  step  by 


F, 


G(c;—u 


RT 


f; 


eq 


(24) 


The  qt  term  in  Equation  23  contributes  to  recovering  the  viscous  heating 
and  compression  work  done  by  the  pressure  and  is  written  as 


Qi 


(Ci-ny 


1 

p 


(-vp+v-n)+(cf-u)-vu 


(25) 


where  the  stress  tensor 


n  =  pu(Vu  +  uV)  (26) 

The  introduction  of  /)  and  gt  seeks  to  repair  an  inconsistent  viscosity 
calculation  among  the  two  distribution  functions,  as  exposed  by  Chapman- 
Enskog  multiscale  expansion  (Peng  et  al.  2003).  However,  because  of  the 
negligible  effect  of  viscous  heating  for  incompressible  flows,  the  following 
simplified  thermal  energy  distribution  proposed  by  Peng  et  al.  (2003)  can 
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be  shown  to  preserve  computational  accuracy  by  eliminating  the  complex 
gradient  calculations  in  Equation  25,  and,  by  extension,  the  necessity  of  ft 
and  gt. 


Peng  et  al.'s  (2003)  simplified  thermal  distribution  functions  follows  the 
same  evolution  as  the  density  distribution  function  and  is  presented  as 

9i  (x  +  cA ,t  +  St]j-gi(x,t)  =—{g,  (x, t)-gr  (x,t))  (27) 

The  new  distribution  function  presented  in  Equation  27  follows  the  same 
collision  and  streaming  evolution  as  the  classic  LBM.  For  the  thermal 

D  RT 

distribution  function,  the  internal  energy  is  represent  by  e  —  — ,  where 
D  =  2  for  two  dimensions  and  is  calculated  as 

P£  =  Yji  (28) 


The  macroscopic  temperature  can  then  be  extracted  by  solving  Equation  28 
for  T.  In  D2Q9,  the  equilibrium  thermal  distribution  is  given  by 


< 


9l!2,3A 


<6,7,8 


2 pe  u2 

~3~V 


pe 

3  3  ci  ■  u 

9(c,.-u) 

3  u‘ 

9 

2  2  c2 

f2  c4 

2  c2 

pe 

3  +  6di'nA 

9(c,.-u)2 

3  u2 

36 

O  T”  O  « 

c 

2  c4 

2  c2 

(29) 


In  the  simplified  thermal  energy  distribution  model,  the  viscosity  is 
consistent  across  both  distribution  function  and  is  given  by 


v  = 


1 

Tv  2 


eft 


(30) 


Another  parameter  of  interest  is  the  thermal  diffusivity,  X,  and  it  is  related 
to  the  thermal  relaxation  time  by 
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X 


(31) 


3.1  Boundary  conditions 

Correct  implementation  of  boundary  conditions  is  especially  important  for 
the  thermal  distribution  function  method  due  to  the  addition  of  a  second 
distribution  function.  In  their  work,  Liu  et  al.  (2010)  applied  a  density 
boundary  condition  shown  along  a  single  wall.  In  this  work,  such  a 
boundary  condition  is  applied  and  is  presented  for  the  case  of  a  general 
wall  as  reformulated  by  the  authors.  In  this  scheme,  the  unknown  density 
distribution  functions  are  given  as 

/„=/;+f  (£■<?)  02) 

where  Q  =  ( Qx  x,  Qy  y)  and  acts  as  a  force  corrector  to  enforce  the  desired 
momentum  on  the  boundary  and  f*  is  a  value  to  be  chosen  and  is  usually 
represented  by  an  approximation  to  fa.  In  this  work,  f*  =  f^q  ,  and  o 
represents  the  set  of  the  three  unknown  distribution  function  indexes  on 
the  wall  such  that  o  =  {a,  /?,  y}.  Following  the  traditional  formulation  for  p 
and  u  =  (it  x,v  y)  and  assuming  that  it  and  v  on  the  wall  are  known,  p,Qx, 
and  Qy  can  be  solved  by  setting  up  the  following  linear  system: 


1  -K  -K 

\ 

P 

if) 

U  -Kc  ~Ky 

QX 

— 

fx 

V  -Ky  ~Ky, 

w 

f 

Jy ) 

f=YJi+YS° 

i  o 

f*  =  YJicix+Yj°c°x 

i  o 

f = y/c  +y^7  c 

Jy  /  1  m  /  ^  °  °y 

i  0 


(33) 


(34) 
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k  =  V w  c 

x  /  j  a  ax 
o 

k  =S^w  c 

y  Z— /  o  oy 

o 

k  =  T w  c2  (35) 

xx  /  j  a  ax  v  7 

o 

k  =  V w  c2 
yy  Z_^ooy 

o 

k  —  w  c  c 

xy  /  j  o  ox  ay 


where  wa  is  the  oth  weight  determined  using  the  D2Q9  template,  fa  -  fa  - 
f*,  and  fa  is  the  yet  to-be-updated  oth  function  value.  It  is  important  to 
note  that  this  formulation  for  fa  on  the  wall  is  derived  using  only  Equation 
32,  p  =  Yiifi  and  pn  =  Si  cji-  Using  Equation  34  and  Equation  35,  values 
for  p,Qx,  and  Qy  can  be  determined  exactly  by  any  preferred  linear  system 

solving  method.  For  the  special  case  of  it  =  v  —  0,  A~x  becomes 


±  KK-KKy  KK-KK 

KXr-K  KKy-K, 


A~1=  0 


0 


yy 


xy 


KKv-K  KK  Ky 


(36) 


xy 


-k 


KX,y  -  Ky  Kkyy  ~Ky 


Note  that  for  all  non-corner  nodes,  kxy  —  0,  thus  yielding  Equation  37 


1 


A-1 


0 


0 


(37) 


Equation  37  is  applied  to  all  walls  in  the  simulation  involving  Rayleigh 
convection.  At  the  corner  nodes,  the  non-equilibrium  bounceback 
condition  must  be  satisfied: 


(38) 
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where  f^eq  =  fteq  -  ft  and  again  ca  =  -cp.  Because  f5  and  f7  for  the  top 
left  and  bottom  right  corners  and  f6  and  /8  for  the  bottom  left  and  top 
right  corners  do  not  stream  into  the  system  (Figure  3),  yet  contribute  to 
the  density  of  the  corner  nodes,  they  are  set  to  their  respective  equilibrium 
function  values. 


Figure  3.  Example  of  non-streaming  lattice 
velocity  directions  in  the  D2Q9  model  for  the 
top  left  channel  corner.  It  follows  that  for  the 
bottom  left  and  top  right  corners,  lattice 
velocities  6  and  8  directions  do  not  stream. 


If  the  temperature  on  the  wall  is  known,  the  thermal  boundary  conditions, 
presented  by  Liu  et  al.  (2010),  are  similar  to  the  density  conditions.  For 
this  case,  the  three  unknown  thermal  distribution  functions  are  given  by 

go=gl  +  woGc  (39) 

The  set  of  distribution  functions,  g*a,  are  determined  arbitrarily  (In  this 
work,  g*a  =  geaq).  Because  the  temperature  on  the  wall  is  known,  it  follows 
that 


ps = RTk = Yji  -  Yj* + Y^g° + w°GJ  (4°) 

1  a  a 


where  Tk  is  the  known  temperature  and  o  again  represents  the  set  of  the 
unknown  distribution  function  indices.  The  quantity  ps*  is  represented  as 
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p£* =Yj9i-Yj9v+Yjl  (41) 

Z  <7  O 


With  Equation  40  K,  Gc  is  determined  by 


G„ 


pe-pe 

E^ 


(42) 


thus  allowing  solutions  of  the  unknown  wall  distribution  functions  from 
Equation  39. 


For  walls  where  the  temperature  is  unknown,  a  second  order  Taylor  series 
expansion  is  applied.  To  obtain  the  temperature  on  the  wall,  the 
temperature  of  the  two  closest  normal  nodes  and  the  local  rate  of  change 
of  the  temperature  in  the  normal  direction  are  required.  On  the  bottom 
wall,  for  example,  the  temperature  is  given  by 


L0  J 


2dT°Js+-z 


3  dy 


Li  j 


2  j 


(43) 


where  8y  is  the  lattice  spacing.  Once  the  temperature  on  the  wall  marker  is 
determined,  the  values  of  the  three  unknown  distribution  functions  can  be 
found  from  Equations  39-42. 

3.2  Validation 

In  this  work,  thermodynamically  driven  convection  was  simulated  to 
demonstrate  computational  validity  for  the  thermal  energy  distribution 
method.  Natural  convection  has  been  shown  to  be  a  conventional  method 
to  simulate  thermal  flows  and  is  an  attractive  validation  case  due  to  ease  of 
implementation  (Nor  Azwadi  et  al.  2006;  Guo  et  al.  2007;  He  et  al.  1998; 
Liu  et  al.  2010;  Peng  et  al.  2003;  Shu  et  al.  1997). 

Here,  a  square  chamber  of  particles,  initially  at  rest,  is  introduced  with  the 
left  wall  maintained  at  a  high  temperature,  T1,  and  the  right  wall 
maintained  at  a  low  temperature,  T0.  The  top  and  bottom  walls  are 
assumed  to  be  adiabatic,  such  that  ^  =  0.  The  temperature  difference 

between  the  two  walls  introduces  a  temperature  gradient  in  the  system. 
Convection  occurs  when  a  newly  introduced  dimensionless  parameter,  the 
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Rayleigh  number,  Ra,  is  above  a  critical  value  causing  the  system  to  tend 
toward  thermal  equilibrium. 

Many  studies  of  Rayleigh  convection  are  examined  under  the  Boussinesq 
approximation,  which  assumes  that  the  thermal  expansion  coefficient,  /?, 
and  the  viscosity,  v,  are  constant.  Under  this  assumption,  the  buoyancy 
force  acting  on  the  system  is  given  as 

pG  =  p(3g0(T-Tm)j  (44) 


where  g0  is  the  acceleration  due  to  gravity,  Tm  is  the  static  average 
temperature  given  by  Tm  =  ,  and/  is  the  unit  vector  pointing  in  the 

opposite  direction  of  gravity.  Two  other  important  quantities  are  the 
dimensionless  Rayleigh,  Ra,  and  Prandtl,  Pr,  numbers  given  as 


Ra 


#7qA  TH3 
vX 


(45) 


Pr  = 


v 

X 


(46) 


where  AT  =  T1  —  T0  and  H  is  the  height  of  the  channel.  To  test  validity  of 
the  method  the  Prandtl  number  is  set  to  Pr  =  0.71  allowing  for  simulation 
of  many  gases.  The  study  is  conducted  over  a  set  of  fixed  Ra.  In  this  work, 
the  incompressibility  constraint  (less  than  5%  mass  loss)  is  determined  by 
the  value  of  / 3g0ATH .  Setting  (3g0ATH  to  a  sufficiently  small  value  will 
maintain  the  system  within  the  incompressible  limit  and  is  chosen  per  the 
Rayleigh  number. 

Once  the  value  of  Pr  and  Ra  is  chosen,  v  and  X  can  be  determined  by 
Equations  45  and  46.  With  v  and  X,  tv  and  rc  are  then  determined  by 
Equations  30  and  31,  respectively. 

The  Nusselt  number,  Nu,  is  another  important  dimensionless  parameter 
that  describes  convective  heat  transfer.  Its  average  value,  Nu,  is  given  by 

0  0 


(47) 
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qx(x,y)  =  uJ(x,y)-X^^l  (48) 

On  the  boundary,  Liu  et  al.  (2010)  velocity  conditions  per  Equation  37  are 
applied  with  all  wall  velocities  set  to  zero.  Similarly,  the  right  and  left 
thermal  boundaries  employ  the  known  temperature  condition  set  to  and 

T0,  respectively.  For  the  top  and  bottom  walls,  an  adiabatic  condition  is 
applied.  The  temperature  is  extrapolated  following  the  second  order 
Taylor  scheme  shown  in  Equation  40.  To  ensure  the  wall  is  adiabatic, 

is  set  to  zero,  where  Tw  is  the  temperature  on  the  walls. 

To  demonstrate  the  validity  of  the  simulation,  the  obtained  results  (Table  2 
are  compared  with  the  simplified  thermal  distribution  function  LBM 
implementation  of  Peng  et  al.  (2003),  and  the  results  from  a  differential 
quadrature  implementation  (noted  as  DQ  in  the  table)  of  the  Navier  Stoke’s 
equations  by  Shu  et  al.  (1997). 


Table  2.  Comparison  of  convection  velocities  for  various  values  of  Ra. 


Ra 

103 

101 x 101 

104 

151 x 151 

105 

201  X  201 

106 

251  x  251 

Umax 

Current 

3.712 

16.275 

33.557 

60.600 

Peng 

3.644 

16.146 

34.261 

63.671 

DQ 

3.649 

16.190 

34.736 

64.775 

y 

Current 

0.810 

0.820 

0.855 

0.855 

Peng 

0.810 

0.820 

0.855 

0.852 

DQ 

0.815 

0.825 

0.855 

0.850 

Vmax 

Current 

3.752 

19.748 

69.058 

224.52 

Peng 

3.691 

19.593 

67.799 

217.57 

DQ 

3.698 

19.638 

68.640 

220.64 

X 

Current 

0.180 

0.120 

0.065 

0.040 

Peng 

0.180 

0.120 

0.065 

0.040 

DQ 

0.180 

0.120 

0.065 

0.035 

Nu 

Current 

1.132 

2.253 

4.524 

8.756 

Peng 

1.117 

2.241 

4.511 

8.737 

DQ 

1.118 

2.245 

4.523 

8.762 
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The  maximum  value  of  the  x  velocity  component,  Umax,  and  its  y  position  at 
a  static  x  location  of  x  =  where  is  L  is  the  length  of  the  channel,  is 
recorded.  In  the  same  way,  the  maximum  value  of  the  y  velocity  component, 
Vmax,  and  its  x  position  at  a  static  y  location  of  y  -  -,  where  H  is  the  height  of 
the  channel,  is  recorded.  The  corresponding  average  Nusselt  number  is  also 
presented  for  each  case.  The  isotherms  and  velocity  streamlines  are 
presented  in  Figure  4  and  Figure  5,  respectively.  The  results  show  very  good 
agreement  with  both  previous  LBM  and  NS  simulations. 


Figure  4.  Isotherms  for  various  values  of  Ra 


Ra=  103 


Ra=  104 
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Figure  5.  Velocity  streamlines  for  various  values  of  Ra. 
Ra=  103  Ra=  104 


Ra=  105 


Ra=  10® 
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4  Microflow  Interactions 

Due  to  the  difficulty  of  simulating  microscale  interactions  with  popular  NS 
based  fluid  solvers  (Lim  et  al.  2002;  Niu  et  al.  2004),  an  extension  to  the 
LBM  able  to  capture  rarefaction  effects  is  introduced  (note  this  implemen¬ 
tation  is  able  to  incorporate  such  effects  solely  by  means  of  an  increase  in 
Kn  and  implementation  of  Diffuse  Scattering  Boundary  Conditions  to 
calculate  the  slip  velocity  with  respect  to  a  solid-fluid  interface  characteristic 
to  microscale  flows).  For  traditional  macroscopic  flows,  rarefaction  effects 
on  the  microscale  are  usually  considered  negligible  and  are  omitted  in 
practice.  To  measure  microscale  interactions,  the  Knudsen  number,  Kn,  is 
employed,  being  a  dimensionless  parameter,  defined  as  the  ratio  of  a 
particle's  mean  free  path  to  the  characteristic  length  of  the  channel.  As  the 
Knudsen  number  increases  (i.e.,  the  mean  free  path  is  on  the  order  of  the 
characteristic  length),  particulate  interactions  must  be  considered.  This 
ensures  that  traditional  Navier-Stokes  solvers  are  generally  unable  to  yield 
accurate  results  in  the  transitional  regime  (0.1  <  Kn  <  1)  due  to  modeling 
a  fluid  as  continuum  with  continuous  macroscopic  quantities.  However, 
microscale  interactions  must  be  studied  with  emphasis  on  rarefaction 
effects  with  discrete  macroscopic  quantities  (Allen  2006).  The  LBM  is  an 
attractive  method  for  microscopic  fluid  flow  study  due  in  part  to  the  lack  of 
an  assumed  continuum  flow  constraint,  as  well  as  a  straightforward 
implementation  of  the  Diffuse  Scattering  Boundary  Conditions,  able  to 
expose  slip  velocity  in  rarefied  gas  flows.  Further,  the  use  of  non-constant 
collision  frequency  allows  direct  implementation  of  microscale  interactions. 

Among  the  rarefaction  effects  to  be  modeled  at  the  microscale  is  the 
existence  of  a  Knudsen  layer  which  exhibits  non-equilibrium  effects  at  the 
boundary  and  spans  about  the  distance  of  one  mean  free  path  from  the 
fluid-solid  interface.  This  layer  is  formed  due  to  the  heightened  frequency 
of  particle  collisions  at  the  boundary  (Homayoon  et  al.  2011)  and  creates  a 
slip  velocity  at  the  fluid-solid  interface  expressed  as  the  particle’s  velocity 
relative  the  boundary  velocity.  For  higher  Kn,  the  effects  of  the  Knudsen 
layer  are  a  dominant  effect  in  the  flow.  Further,  collisions  at  boundaries 
are  reflected  in  a  manner  slightly  inconsistent  with  the  equilibrium  state 
described  by  the  Maxwellian  distribution  function.  Certain  treatments 
must  be  implemented  for  accurate  simulation  of  a  Knudsen  layer,  as  well 
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as  boundary  conditions  able  to  describe  fluid-boundary  collisions.  (Lilley 
and  Sader,  2008) 

A  variety  of  proposals  exist  to  effectively  tie  the  LBM  with  accurate 
microchannel  representation.  Among  these  are  the  introduction  of  effective 
quantities  as  employed  by  Ghazanfarian  and  Abbassi  (2010)  and  Homayoon 
et  al.  (2011),  and  are  written  in  terms  of  the  originally  calculated  quantities. 
Hoyamoon  et  al.  (2011)  proposed  an  effective  dynamic  viscosity  that  yields 
an  effective  relaxation  time  where  both  parameters  were  inversely 
proportional  to  some  function  of  Kn.  Lim  et  al.  (2002)  introduced  an 
updated  Kn  inversely  proportional  to  the  pressure  distribution  in  the 
channel  effectively  updating  the  relaxation  time  due  to  the  assumption  that 
the  channel  pressure  is  proportional  to  the  mean  free  path  and  the  mean 
free  path  is  proportional  to  the  relaxation  time.  Niu  et  al.  (2004)  presented 
tv  oc  Kn  modified  by  the  specific  heat  ratio  for  a  monatomic  gas. 

Of  greatest  significance  though,  is  that  most  authors  agree  that  the 
dimensionless  relaxation  time  is  proportional  to  the  Knudsen  number 
since  tv  oc  Kn  (Ghazanfarian  and  Abbassi  2010;  Homayoon  et  al.  2011; 

Lim  et  al.  2002;  Niu  et  al.  2004).  Therefore,  the  introduction  of  a  fixed  Kn 
is  directly  reflected  in  an  update  of  the  collision  frequency. 

Letting  X  represent  the  mean  free  path  of  the  considered  particles  and  H 
be  the  characteristic  length  of  the  channel,  here  representing  the  channel 
height,  the  Knudsen  number  is  given  by 

Kn =A  (49) 


The  mean  free  path  can  then  by  calculated  by  X  =  ( v)Q ,  where  (v)  is  the 
average  magnitude  of  velocity  of  the  particles  enclosed  in  the  system  given 

1 8  RT 

by  (v)  =  —  and  0  is  the  relaxation  time  related  to  xv  by  0  =  r vSt.  (Niu  et 

al.  2004)  In  the  D2Q9  model,  ( v )  is  transformed  by  c  —  V3 RT  to  (v)  — 


c_  — .  Thus,  writing  X  in  terms  of  rv  yields 


(50) 
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The  characteristic  length  is  then  given  by 

H  =  Nh8x  (51) 

where  NH  is  the  number  of  lattice  sites  that  span  the  width  of  the  channel 
and  8X  is  the  lattice  site  spacing  and  is  assumed  to  be  constant  throughout. 
From  Equation  50  and  Equation  51,  the  evaluation  of  Kn  can  be  updated, 
yielding: 


Kn  = 


,  IjE  tv6' 

v  3.1  Nh6x 


(52) 


Noting  that  c  —  -f- and  implementing  a  correction  factor  of  —  -  to  maintain 

ot  2 

second  order  accuracy  (in  accordance  with  with  the  Chapman-Enskog 
analysis  (Ghazanfarian  and  Abbassi  2010)),  Equation  52  becomes 


Kn 


(53) 


The  Kn  has  now  been  effectively  tied  to  relaxation  time  as  assumed  by  the 
present  authors.  With  this  formalization  of  Kn  for  the  D2Q9,  the 
relaxation  time  can  now  be  determined  in  terms  of  a  desired  Kn.  With  the 
calculated  r,  A  can  also  be  computed  by  Equation  50.  In  order  to  ensure  a 
more  accurate  Knudsen  layer,  Ghazanfarian  et  al.  (2010)  proposed  a  wall 
function  that  modifies  A  in  Equation  50  by 


A 


A 


1  +  ip 


+  W 


H-y 


A 


(54) 


where  ip  is  given  by  ip(x)  —  e~Cx,  C  is  a  constant  set  to  unity  in  this  work, 
P/’  is  the  previously  introduced  Prandtl  number  and  y  is  the  normal 
distance  relative  to  the  bottom  wall.  Further,  Ghazanfarian  and  Abbassi 
(2010)  also  formulated  a  modification  to  tv  in  attempt  to  couple  the 
contributions  generated  by  the  density  and  internal  energy  distribution 
functions  that  is  given  by 


ERDC  TR-14-6;  Report  2 


25 


T'  ~2  Pref  | 

( 

6 

1  +  IP 

(  \ 

y 

UJ 

+  w 

\H~y ) 
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P  \ 

.  €ref  , 

1 

+  2 


(55) 


Note  that  t'v  is  written  in  terms  of  tv,  e,  and  p.  This  calculation  considers 
the  desired  microscale  interactions  through  the  implementation  of  tv 
derived  from  a  chosen  Kn,  internal  energy  updated  by  thermal 
interactions,  and  pressure  contributions  due  to  the  linear  relationship  of 
pressure  and  density  in  the  LBM.  Implementation  of  Equation  55  couples 
density  and  internal  energy  interactions  in  LBM  allowing  for  both  and 
pressure  and  thermal  interactions  to  be  simulated  on  the  microscale. 

Since  Pr  —  p  v  oc  tv,  and  X  oc  rc,  the  thermal  relaxation  time  can  then  be 
computed  in  terms  of  Pr  as 


(56) 


4.1  Boundary  conditions 

To  correctly  model  the  slip  velocity  on  the  walls,  the  standard  bounceback 
condition  is  shifted  to  a  discretized  Maxwell  boundary  condition,  referred 
to  as  the  DSBC.  The  DSBC  as  presented  (Tang  et  al.  2005)  is  shown  and 
was  implemented  in  this  work. 

With  knowledge  of  the  wall  velocity,  the  motion  of  the  particles  near  the 
wall  can  be  expressed  with  slip  reference  velocity  of  =  ct  -  uw.  The 
reference  velocity  projection  over  the  unit  normal  to  the  wall  is  then  given 
as  ■  n,  where  n  is  the  normal  unit  vector  pointing  into  the  flow.  With  the 
values  of  ft  ■  n,  the  respective  ft  can  be  grouped  into  three  sets.  Particles 
with  values  of  £  ■  n  <  0  must  be  incident  to  the  wall  and  are  grouped  in  / 
Particles  with  &  •  n  >  0  are  reflected  particles  and  are  placed  in  Ir.  Finally, 
particles  with  &  ■  n  =  0  must  be  "grazing"  particles  and  are  placed  in  Ig. 

In  this  grouping,  information  about  incident  and  grazing  particles  are 
already  known  and  the  reflected  particles  must  be  solved  for.  According  to 
kinetic  theory,  the  reflected  particle  distribution  function  is  related  to  that 
of  the  incident  particles  by  (Struchtrup  2013) 
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fi, 


jeir 


7 Z(&—>  c, 


Zr* 


(57) 


Here  31  (q  ->  c')  is  the  scattering  kernel  which  represents  the  probability 
that  particles  which  hit  the  boundary  with  velocity  { c k,  c'k  +  dc '}  will  return 
to  the  flow  with  velocity  { ck ,  ck  +  dc'}.  Due  to  the  consideration  of  only 
nonporous  and  nonabsorbent  boundaries,  the  normalization  constraint  for 
3l(ci  ->  c-)  yields  (Struchtrup  2013;  Tang  et  al.  2005) 


hot 

J  Tl(ci  — ■>  c^dc  =  1 

ielt 


(58) 


which  guarantees  that  all  boundary  colliding  particles  are  returned  to  the 
flow.  For  purposes  in  the  LBM,  discretization  of  Equation  57  is  written  as 

(C  ■*)/:  <59> 

ielt 


In  attempt  to  represent  multiple  particle-boundary  interactions,  various 
formulations  of  31  are  considered  by  kinetic  theory.  In  his  model,  Maxwell 
chose  to  consider  the  specular  and  diffusive  derivations  of  31  and  couple 
them  linearly  with  the  tangential  momentum  accommodation  coefficient, 
ov,  such  that  ov  e  [0, 1]  by  applying  ov  to  the  diffuse  condition  and  ov-  1 
to  the  specular  condition.  Considering  the  specular  and  diffusive  boundary 
conditions,  the  DSBC  in  the  LBM  may  be  expressed  as 


<7, 


fj  ( * ,  t T  St )  =  ( 1  —  a  v )  ff  ( x ,  t  +  St )  ^ j  •  n )  + 

;fjq(x,PwA»d  +  8t 


E 

i£li  ( 

)/,\ 

[x,f  +  <5f) 

)/r( 

x,t 

jw,uw,t  +  8tj 

(60) 


where  j  E  Ir.  Note  that  taking  ov  =  1  considers  the  particle  reflection  to  be 
completely  diffusive  and  ov  =  0  takes  the  particle  reflection  to  be  completely 
specular.  In  this  work,  the  value  of  ov  —  1  is  exclusively  considered.  The 
diffusive  interaction  ensures  that  incoming  particles  are  reflected  in  such  a 
way  that  the  particle  "loses  its  memory"  of  its  incoming  direction  and  is 
after  scattered  such  that  both  the  tangential  and  normal  components  of  its 
direction  are  altered.  (Niu  et  al.  2006;  Struchtrup  2013;  Tang  et  al.  2005)  It 
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has  also  been  shown  by  Niu  et  al.  (2006)  that  for  steady  flow  where  the 
velocity  varies  only  the  streamwise-normal  direction  the  slip  velocity  may 
written  in  terms  of  only  Kn  and  u. 

4.2  Validation 

Correct  implementation  of  microchannel  support  for  the  LBM  was 
validated  through  simulation  of  Poiseuille  flow  with  heated  walls  over  a  set 
of  fixed  Kn.  A  strong  pressure  gradient  of  ratio  2  across  the  channel  is 
implemented  using  the  enforced  density  boundary  condition  proposed  by 

p 

Zou  et  al.  (1997).  In  the  work,  —  =  2,  with  P1  being  the  inlet  pressure  and 

A 

P2  is  the  outlet  pressure.  The  velocity  is  then  calculated  in  order  to 
maintain  the  desired  pressure.  At  the  north  and  south  walls,  the  velocity 
condition  in  Equation  37  was  first  applied  to  ensure  the  wall  velocity  was 
zero  and  then  DSBC  conditions  were  applied.  Due  to  the  DSBC,  the 
presence  of  a  slip  velocity  was  observed  as  expected  for  microscale  flows. 

The  bottom  and  top  walls  have  an  imposed  internal  energy  of  and  e2, 

respectively,  with  an  internal  energy  ratio  corresponding  to  —  =  2 .  At  the 

e2 

inlet,  the  internal  energy  is  as  measured  with  respect  to  the  top  wall  is  —  = 

e2 

1.5.  To  enforce  the  internal  energy  on  the  walls,  the  thermal  boundary 
condition  presented  in  Equations  39-42  are  used  by  first  determining  the 
known  temperature  from  the  internal  energy  with  e  —  RT. 

To  enforce  the  effects  of  the  Knudsen  layer,  the  wall  function  modification 
of  A  shown  in  Equation  54  is  implemented.  The  adjustment  of  tv  by 
Equation  55  is  also  implemented,  effectively  coupling  the  microchannel 
and  thermal  interactions.  With  the  use  of  Equation  55,  the  viscosity  and 
thermal  diffusivity  are  now  both  dependent  on  pressure  and  internal 
energy.  Due  to  this,  it  can  be  observed  that  the  pressure  distribution  is  no 

longer  linear.  The  pressure  is  presented  after  being  normalized  by - - 

P2 

where  Pt  is  the  expected  linear  pressure  distribution  given  in  terms  of  x*  = 
j  with  L  being  the  length  of  the  channel  by  Pt  =  P1+  x*(P2  -  Pj.  Deviation 

from  linear  pressure  is  compared  to  the  results  found  by  Ghazanfarin  and 
Abbassi  (2010)  and  is  presented  for  various  Kn  in  Figure  6. 

To  demonstrate  the  validity  of  the  implemented  model,  comparisons  of 
velocity  profiles  normalized  by  the  bulk  velocity  are  presented  for  various 
Kn.  The  presented  results  are  compared  with  those  produced  by 
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Ghazanfarian  and  Abbassi  (2010)  and  are  shown  in  Figures  6-9.  Due  to 
thermal-microscale  coupling,  the  velocity  is  updated  depending  on 
pressure,  internal  energy,  and  Knudsen  layer  modeling,  making  the  velocity 
profiles  an  interesting  calculation.  It  is  seen  that  the  yielded  results  agree 
very  well  with  those  presented  by  Ghazanfarian  and  Abbassi  (2010). 


Figure  6.  Deviation  from  linear  pressure  distribution  for  various  Kn. 


Current  Kn  =  0.1 


Ghazanfarian  et 
al.  Kn  =  0.1 


Current  Kn  =  0.2 


Ghazanfarian  et 
al.  Kn  =  0.2 


Current  Kn  =  0.4 


Ghazanfarian  et 
al.  Kn  =  0.4 


Figure  7.  Normalized  velocity  profile  corresponding  to  Kn=  0.1. 
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Figure  8.  Normalized  velocity  profile  corresponding  to  Kn=  0.2. 


Figure  9.  Normalized  velocity  profile  corresponding  to  Kn=  0.4. 
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5  Summary 

This  report  demonstrated  the  validity  of  various  physical  extensions 
implemented  within  the  LBM.  It  was  shown,  through  a  variety  of  worked 
examples,  that  the  LBM  is  able  to  accurately  account  for  flows  around 
complex  immersed  boundaries,  thermal  convection  driven  flows,  and 
microchannel  flows.  Further,  this  effort  also  served  to  set  the  LBM  apart  as 
a  computationally  accurate  method,  especially  for  microscale  flow 
scenarios.  In  particular,  this  accuracy  was  demonstrated  in  simulations  of 
(l)  flow  around  an  immersed  cylinder,  (2)  natural  convection  in  a  square 
channel,  and  (3)  pressure  driven  flow  with  heated  walls  and  Knudsen  layer 
accommodation.  Of  particular  interest  is  the  effective  coupling  of  thermal 
and  microchannel  support.  Given  the  present  work  as  a  foundation,  future 
efforts  may  focus  on  particulate  flows  through  microscale  electronics, 
deformation  and  heat  transfer  of  grain-structured  materials,  and 
rarefaction  effects  on  microscale  boundaries. 
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Appendix  A:  Discretized  Delta  Function 
Validation 

Consider  a  Lagrange  marker  located  at  Xs  -  (xj  +  Ax,  y;-  +  Ay)  where  ftj  = 
(xj,  yj)  is  the  nearest  lower  left  Eulerian  marker  relative  to  Xs  and  xt  > 

Ax  >  xi+1,  yj  >  A y  >  yj+1.  If  the  force  at  the  Eulerian  nodes  is  updated  by 

Equation  8  and  applied  to  a  domain  H,  the  magnitude  of  force  distributed 
to  rLj  is  given  by 


YPh  (i  -  n  + 1  +  Ax)  ^  (j~n  +  l  +  Ay)  (Al) 

i= 1  j= 1 

Let  a  be  the  first  summation  and  /?  the  second  summation,  such  that  a  — 
Yi=1 8h(i-n  +  1  +  Ax)  and  /?  =  jj=1 8h(j  -  n  +  1  +  Ay).  Note  that  terms 

above  n  =  4  equal  o.  Expanding  a,  we  find 

Q  =  6)  (— 2  +  Ax)  +  8h  (— 1  +  Ax )  +  Sh  (Ax)  +  8h  (l  +  Ax )  (A2) 


Evaluating  first  term  of  a  yields 


«„(-2  +  Ax)=- 


1  +  cos 


[-271 

Axji  ' 

I 

2 

r 

2  ) 

(A3) 


Applying  the  trigonometric  identity  for  cos  (a  +  b),  we  find  Equation  A3 
evolves  into  Equation  A4 


S„(- 2  +  A*)  =  i 


,  —2 71  A xn  .  —2 n  .  A xn 

1  +  cos - cos - sin - sin - 


(A4) 


After  following  this  procedure  of  evaluating  the  delta  function  and 
applying  the  trigonometric  identity  for  cos  (a  +  b  )  for  the  remaining  three 
terms  in  a,  we  group  terms  and  find 
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4  + cos 


A xji 


sin- 


2 
A xji 


—  2  71  — 71  On  7l' 

cos - b  cos - b  cos - b  COS  — 

[  2  2  2  2) 

.  —2 n  .  — 7i  .  On  .  n' 

sin - b  sin - b  sin - b  sin — 

2  2  2  2 


(A5) 


Evaluating  the  four  inner  cosines  and  sines  proves  that  a  =  l  for  all  Ax. 
This  same  procedure  can  be  applied  for  /?  and  it  becomes  apparent  that  /? 
=  l  as  well.  Therefore,  a/3  =  l  for  enclosed  discretized  Eulerian  markers 
independent  of  the  location  of  the  forcing  Lagrange  marker  at  Xs. 
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Appendix  C:  Definition  of  Mathematical 
Symbols  Used  in  IMB 


Table  Cl.  Definition  of  Mathematical  Symbols  used  in  IMB. 


Symbol 

Definition 

Units 

drag  coefficient 

Fr 

force  reapplied  to  flow 

Fs 

immersed  boundary  body  force 

h 

immersed  boundary  moment  of  inertia 

T0 

immersed  boundary  torque 

X 

Lagrangian  node  position 

m 

mass 

V 

^  cm 

object  center  of  mass 

Du 

2D  lattice  force  mollifier 

Greek 

(A) 

angular  velocity 

Sh 

force  distribution  mollifier 
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Appendix  D:  Definition  of  Mathematical 
Symbols  Used  in  Thermal  Interactions 


Table  Dl.  Definition  of  Mathematical  Symbols  used  in  Thermal  Interactions. 


Symbol 

Definition 

Units 

do 

acceleration  due  to  gravity 

T 

1  m 

average  temperature 

R 

gas  constant 

T 

temperature 

9eq 

thermal  equilibrium  distribution  function 

q 

viscous  heating  and  compression  contribution 

AT 

wall  temperature  difference 

Greek 

£ 

internal  energy 

n 

stress  tensor 

X 

thermal  diffusivity 

unitless 

Tc 

thermal  distribution  relaxation  time 

unitless 

P 

thermal  expansion  coefficient 

unitless 
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Appendix  E:  Definition  of  Mathematical 
Symbols  Used  in  Microflow  Interactions 


Table  El.  Definition  of  Mathematical  Symbols  used  in  Microflow  Interactions 


Symbol 

Definition 

Units 

Nh 

characteristic  length 

31 

scattering  kernel 

Greek 

? 

fluid-wall  relative  velocity 

X 

mean  free  path 

unitless 

0 

relaxation  time 

s1 

tangential  momentum  accommodation  coefficient 

unitless 
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